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Background: Dysfunctional integration of distributed brain networks is believed to be the 
cause of schizophrenia, and resting-state functional connectivity analyses of schizophrenia 
have attracted considerable attention in recent years. Unfortunately, existing functional 
connectivity analyses of schizophrenia have been mostly limited to linear associations. 

Objective: The objective of the present study is to evaluate the discriminative power of 
non-linear functional connectivity and identify its changes in schizophrenia. 

Method: A novel measure utilizing the extended maximal information coefficient was 
introduced to construct non-linear functional connectivity. In conjunction with multivariate 
pattern analysis, the new functional connectivity successfully discriminated schizophrenic 
patients from healthy controls with relative higher accuracy rate than the linear measure. 

Result: We found that the strength of the identified non-linear functional connections 
involved in the classification increased in patients with schizophrenia, which was opposed 
to its linear counterpart. Further functional network analysis revealed that the changes of 
the non-linear and linear connectivity have similar but not completely the same spatial 
distribution in human brain. 

Conclusion: The classification results suggest that the non-linear functional connectivity 
provided useful discriminative power in diagnosis of schizophrenia, and the inverse but 
similar spatial distributed changes between the non-linear and linear measure may indicate 
the underlying compensatory mechanism and the complex neuronal synchronization 
underlying the symptom of schizophrenia. 



Keywords: schizophrenia, resting-state functional connectivity, non-linear, extended maximal information 
coefficient, compensatory 



INTRODUCTION 

Schizophrenia, which is characterized by delusions, auditory hal- 
lucinations, and impairments in memory, attention, and execu- 
tive function, is one of the most devastating, cryptic, and costly 
psychiatric disorders (van Os et al, 2010; Sui et al., 2012). 
It brings not only great suffering to patients but also signifi- 
cant costs to society. Traditionally, the diagnosis of schizophre- 
nia depends on the observation of psychiatric symptoms and 
longitudinal courses, while modern diagnoses of such psychi- 
atric disorders require objective neurological measures (Kendler, 
2009; Insel et al, 2010; Shen et al, 2010). The underlying eti- 
ology and mechanisms of schizophrenia are still unclear, but 
it is believed that the dysfunctional integration of distributed 
brain networks leads to this mental disease (Friston and Frith, 
1995; Andreasen et al., 1998). Thus, based on fMRI data, func- 
tional connectivity research on schizophrenia has attracted con- 
siderable attention in recent years (Camchong et al, 2011; 
Pettersson-Yeo et al., 2011; Bassett et al., 2012; Fornito et al., 
2012; Liu et al, 2012; Tu et al, 2012; Vertes et al., 2012; 
Zalesky et al, 2012). Furthermore, the introduction of multivari- 
ate pattern classification to behavioral and cognitive neuroscience 



increases the potential of functional connectivity in clinical diag- 
noses of this mental disease (Shen et al., 2010; Fan et al., 
2011). 

Functional connectivity, defined as statistical associations 
between remote neurophysiological events, aims to character- 
ize the communication between different brain regions (Friston 
et al, 1993). Most functional connectivity analyses use temporal 
correlations or covariance to examine the simultaneous coupling 
between two time series. Fluctuations in the blood oxygen level 
dependence (BOLD) signal have attracted attention since the 
1990's (Biswal et al, 1995, 1996). In contrast to functional con- 
nectivity, which need previous acknowledge about the disease and 
did not estimate the potential functional connectivity of other 
ROIs (Lowe et al, 1998, 2000; Xiong et al, 1999; Cordes et al, 
2000; Hampson et al, 2002), whole brain functional connectivity 
analyses of schizophrenia which require less field specific knowl- 
edge have attracted considerable attention in recent years (Liang 
et al, 2006; Liu et al, 2006; Lynall et al, 2010; Alexander-Bloch 
etal, 2013). 

Previous whole brain functional connectivity analyses of 
schizophrenia which primarily used temporal correlation or 
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covariance as the measure limited to consideration of only the 
linear associations (Liang et al., 2006; Liu et al, 2006; Lynall 
et al., 2010; Alexander-Bloch et al, 2013). Because correlation 
calculations in the full-lag space are computationally expensive 
(Cecchi et al., 2007), correlation related studies have prevalently 
used the zerotfi lag correlation. The hemodynamic response can 
significantly reduce the computational expense for its limited 
duration (Li et al., 2009), but presents another problem, namely, 
the hemodynamic response function (HRF) varies across sub- 
jects and brain regions (Buckner et al, 1998; Miezin et al., 2000; 
Lee et al, 2001; Saad et al, 2001), and the zerof/i lag correla- 
tion is sensitive to changes in regional HRFs. Specifically, the 
simple correlation (zerof/i lag) is sensitive to the lag between 
time series that cannot sufficiently depict the functional inter- 
actions of the human brain (Smith et al., 2011). Thus, we 
believe that the underlying neural activity cannot be accu- 
rately reflected by the zerot/i lag correlation- or covariance-based 
methods. 

Investigating functional connectivity in schizophrenia from 
the frequency domain is an alternative of the temporal corre- 
lation (Fallani et al, 2010; Salvador et al, 2010). Although the 
spectral analogs of functional connectivity such as coherence 
fixed some of the problems (i.e., lag between time series) faced 
by the simple correlation methods, coherence explores only lin- 
ear relationships between time series (Sun et al., 2004; Smith 
et al., 2011). However, the non-linearity of the HRF has been 
reported by several studies (Buxton et al, 1998; Friston et al., 
2003; de Zwart et al, 2009; Daunizeau et al, 2012). As an out- 
put of HRF, the BOLD signal also has non-linear properties 
(Vazquez and Noll, 1998; Xie et al, 2008a,b; Zhang et al, 2008). 
Notably, non-linear relationships between time series extracted 
from resting state BOLD signals have also been confirmed 
(Lahaye et al., 2003). 

The first goal of this study is to evaluate the dis- 
criminative power of non-linear functional connectivity in 
schizophrenia which may have potential application in diag- 
nosis of neuropsychiatric disorders. The second goal is to 
investigate the changes of the non-linear functional connec- 
tivity in schizophrenia which may contribute to the patho- 
physiology of schizophrenia. Aiming at exploring non-linear 
associations in schizophrenia, we present a novel method 
that introduces the extended maximal information coefficient 
(eMIC) to construct functional connectivity of the whole 
brain. 

The eMIC represents the difference between the MIC which 
is a statistical method for detecting various associations between 
pairs of variables in large data sets (Reshef et al., 2011) and 
the square of the Pearson correlation coefficient (PCC). In a 
previous study, the square of PCC was proved to be equal 
to the Hilbert-Schmidt Independence Criterion (HSIC) and 
gave excellent performance in feature selection (Song et al., 
2007). By applying the MIC, PCC and eMIC with fMRI 
data, we investigated the discriminative power of non-linear 
and linear functional connectivity. Based on the classifica- 
tion result, we further evaluated the changes and spatial dis- 
tribution of linear and non-linear functional connectivity in 
schizophrenia. 



MATERIALS AND METHODS 
MATERIALS 

Participants 

The fMRI data used in this study were acquired from 64 
participants who were all right-handed native Chinese speakers. 
Participants consisted of 32 patients suffering from schizophre- 
nia and 32 healthy controls. All the patients were recruited from 
outpatient departments and inpatient units at the Department of 
Psychiatry, Second Xiangya Hospital of Central South University 
in Changsha, China, between March 2006 and October 2007, 
and satisfied the Diagnostic and Statistical Manual of Mental 
Disorders, Fourth Edition (DSM-IV; American Psychiatric 
Association, 1994) criteria. Five patients were medication-free, 
while the others accepted antipsychotic drugs at the time of 
image acquisition. The healthy subjects were recruited by adver- 
tisements, and matched to the patients on age and gender. 
None of them had major head trauma, history of alcohol or 
drug dependence, or history of neurological disorder. Written 
informed consents were obtained from all the subjects who took 
part in this study. This study was approved by the Medical 
Research Ethics Committee of the Second Xiangya Hospital, 
Central South University. Details about the participants were 
displayed in Table 1 . 

Imaging protocol 

In the experiments, subjects were instructed simply to keep their 
eyes closed, to relax, remain awake and perform no specific cog- 
nitive exercise. After each session, subjects were asked whether 
they were awake in the previous session and all the subjects con- 
firmed. MRI scans were performed with a 1.5T GE Signa System 
(GE Signa, Milwaukee, Wisconsin, USA). To reduce the head 
movements, subjects' heads were fixed by using foam pads with 
a standard birdcage head coil. The functional MRI images were 
collected by using a gradient-echo echo planar imaging sequence. 
The imaging parameters were as follows: repetition time/echo 
time = 2000/40 ms, thickness/gap = 5/1 mm, field of view = 
240 x 240 cm, flip angle = 90°, matrix = 64 x 64, and slices = 
20. Each functional resting state session lasted ~6min, and 180 
volumes were obtained. 

Data preprocessing 

Data preprocessing was performed by using statistical paramet- 
ric mapping software package (SPM2, Wellcome Department 
of Cognitive Neurology, Institute of Neurology, London, UK, 
http://www.fil.ion.ucl.ac.uk/spm/). In each subject, the first 5 



Table 1 | Characteristics of participants in this study. 



Variables 


Mean + SD 


P-value 




Patient 


Control 




Sample size 


32 


32 




GenderflWF) 


25/7 


23/9 


0.85 


Age(years) 


24 ± 5.66 


25.01 ±4.50 


0.92 


PANSS 


80.06 ± 16.55 







PANSS: Positive and Negative Syndrome Scale 
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volumes of the scanned data were discarded for magnetic sat- 
uration effects. The remaining 175 volumes were corrected by 
registering and reslicing for head motion. Next, the volumes 
were normalized to the standard echo planar imaging template in 
the Montreal Neurological Institute (MNI) space. The resulting 
images were spatially smoothed with a Gaussian filter of 8 mm 
full-width half-maximum kernel, detrended to abandon linear 
trend and then temporally filtered with a Chebyshev band-pass 
filter (0.01-0.08 Hz). The registered fMRI volumes with the MNI 
template were further divided into 116 regions according to the 
automatic anatomical labeling atlas (Schmahmann et al., 1999; 
Tzourio-Mazoyer et al, 2002). The atlas divides the cerebrum into 
90 regions (45 in each hemisphere) and divides the cerebellum 
into 26 regions (nine in each cerebellar hemisphere and eight in 
the vermis). All region of interest masks were generated using 
the free software WFU_PickAtlas (version 3.0, http://www.ansir. 
wfubmc.edu) (Maldjian et al, 2003). 

Regional mean time series were obtained for each individual by 
averaging the functional MRI time series over all voxels in each of 
the 116 regions. Each regional mean time series was further cor- 
rected for the effects of head movement by regression on the time 
series of translations and rotations of the head estimated in the 
course of initial movement correction by image realignment. The 
residuals of these regressions constituted the set of regional mean 
time series used for functional connectivity analysis (Achard et al., 
2006). 

We evaluated functional connectivity between each pair of 
regions using the MIC, PCC, and eMIC. Thus, for each subject, 
we obtained three resting state functional networks captured by 
three 116 x 116 symmetric matrixes respectively. Removing 116 
diagonal elements, we extracted the upper triangle elements of 
the functional connectivity matrix as classification features, i.e., 
the feature space for classification was spanned by the (116 x 
1 15)/2 = 6670 dimensional feature vectors. 

METHODS 
MIC and emic 

In this section, we provide a brief description of the MIC and 
eMIC for detecting the association between two time series. Two 
time series can be viewed as a set of ordered data pairs. The MIC 
between a set of ordered pairs will not change if the rank of pairs 
is disturbed but the relative ranks of the x- and y-values do not 
change. If two variables are related to each other, then a grid can 
be drawn on the scatterplot of the two variables that encapsulates 
that relationship. Based on this concept, this method investigated 
all the grids up to the maximal grid resolution, which depends on 
the sample size (Reshef et al., 2011). The MI is defined as follows: 

MI(X; Y) = H(X) - H(X\ Y) ( 1 ) 

= H(Y) - H(Y\X) 
= H(X) + H(Y) - H(X, Y) 

where H(X) and H(Y) are the marginal entropies, H(Y \X) and 
H(Y \X) are the conditional entropies, and H(X, Y) is the joint 
entropy of X and Y. 



Formally, let G xxy be all the possible partitions with rows x 
and y columns (width of rows are different; width of the columns 
are different, too) of the scatterplot for the ordered pairs of two 
vectors. As the partitions were not equal, there are many possible 
partitions with x rows and columns of the scatterplot, let I g denote 
the MI for one possible partition with x x y grids that are applied 
to the ordered samples of the two vectors. For fair comparison 
between grids of different resolution, the MI values of different 
partitions with x x y grids of scatterplot for the ordered pairs of 
two vectors are normalized to the interval [0,1]. The m xxy is 
defined as 

m xxy = max {7^}/logmin{x,y} (2) 

MIC is the maximal value of m xxy over the ordered pairs (x, y), 
x < n, y < n, n is the length of the vectors (i.e. number of fac- 
tors in the vector). In practice, to accelerate computation, it is 
not necessary to compute m xxy over all (x,y), x = n, y = n. 
Alternatively, we can compute the m xxy over all (x,y) xy < B. 
Empirically, B is defined as B = n 0 6 . We can then define 

MIC = max{m xx y ] (3) 

xy<B 

Additionally, the eMIC can be defined as 

eMIC = MIC - HSIC = MIC - PCC 2 (4) 

Refer to (Song et al., 2007; Reshef et al., 201 1) for more details of 
the MIC and HSIC. 

Identification of Features with High Discriminative Power 

Due to the noise, registration error, and inter- individual anatomi- 
cal differences, only a small number of the 6670 features are highly 
discriminative. To achieve good discriminative performance as 
well as resistance to noise and individual disparity, the first step 
of constructing the classification model was selecting those fea- 
tures with high discriminative power to construct the feature 
space for classification. The discriminative power of a feature 
can be quantitatively measured by its relevance to classification 
(Guyon and Elisseeff, 2003). Here, we used the Kendall tau rank 
correlation coefficient (Kendall and Jean, 1990) which provides 
a distribution free test of independence between two variables. 
Let (xi,yi) , (x2, yi) , . ■ . , (xn, yti) be a set of samples of the 
joint random variables X and Y respectively. Any pair of samples 
(Xi, yn and (xj, yj) are said to be concordant if the ranks for both 
elements agree. The Kendall tau coefficient is defined as: 

(number of concordant sample pairs) — 

(number of discordant sample pairs) 

x = (5) 

number of total sample pairs 

For any sample (%;, yA , i = 1, 2, . . . , N of the variables need only 
consider the pairs between itself and the rest N — 1 samples. 
When i was changed from 1 to N, each pairs was counted twice. 
Then, number of total simple pairs = 0.5 x N x (N — 1). 

Suppose there are m samples in the control group and n 
samples in the patient group. Let xy denote the functional con- 
nectivity feature i of the ;'th samples and y; denote the class label 
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of this sample (+1 for controls and —1 for patients). The Kendall 
tau correlation coefficient of the functional connectivity feature i 
can be defined as 



m x « 

Because the samples in the same groups are neither con- 
cordant pairs nor discordant pairs, the relationship between 
two samples that belong to the same group is not consid- 
ered, number of total simple pairs = (N — m) x (N — n) , N = 
m + n. The n c and are the numbers of concordant and dis- 
cordant pairs between the two groups, respectively. For a pair of 
two-observation data sets [Xij,yj] and {xj k ,yk}, it is neither con- 
cordant pair nor discordant pair if y ( - = y k ; it is a concordant pair 
when 

sgn(x, ; - - x ik ) = sgn(y, - y k ) (7) 

where sgn( ) is a signum function. Correspondingly, it is a discor- 
dant pair when 

sgn(xy - Xik) = -sgn(>/ - yk) (8) 

Thus, a positive correlation coefficient t, represents that the ith 
functional connectivity exhibits significant decrease in the patient 
group compared to the control group, while a negative x, repre- 
sents that the z'th functional connectivity increases in the patient 
group. The discriminative power was defined as the absolute value 
of the Kendall tau correlation coefficient. Then we ranked all 
the features according to their discriminative power and selected 
those with correlation coefficient over a threshold as the final 
feature set for classification. 

Since a leave-one-out cross-validation strategy was introduced 
to estimate the generalization ability of the classifiers (see below) 
and the training data set for feature ranking is slightly different 
in each iteration of the cross-validation, the first selected fea- 
tures differed slightly from iteration to iteration. Therefore, the 
contribution of different regions to classification was not evenly 
distributed, and some regions formed many highly discrimi- 
nating functional connections with other regions, while some 
did not. Consensus functional connectivity was introduced here, 
which was defined as the functional connectivity feature appear- 
ing in each cross-validation iteration. Region weight, representing 
the relative contribution to the identification of schizophrenic 
patients, was denoted by its occurrence number in the consensus 
functional connections in this study. The consensus functional 
connectivity discriminative power was denoted by the mean of its 
discriminative powers across all iterations of the cross-validation. 

Support vector classification and performance evaluation 

When the data set of features with high discriminative power 
was obtained, support vector machines (SVM) with linear ker- 
nel function were employed to solve the classification problem 
(Vapnik, 1995; Bishop, 2006). Due to the classification results 
were influenced by the number of involved features, classification 
accuracies with fixed parameter setting of SVM (LIBSVM3.il: 
Linear kernel, C = 1 ) using a wide range of feature number were 
reported. Due to our limited number of samples, we used a leave- 
one-out cross-validation strategy to estimate the generalization 



ability of our classifier. The performance of a classifier can be 
quantified using the generalization rate (GR, i.e., the rate of 
all the subjects correctly classified), sensitivity (SS, i.e., the rate 
of the patients correctly classified) and specificity (SC, i.e., the 
rate of the controls correctly classified) based on the results of 
cross-validation. 

RESULTS 

We constructed whole brain functional connectivity using 
the PCC, MIC, and eMIC based on fMRI data collected 
from schizophrenic patients and matched healthy controls. 
Multivariate pattern classification was then introduced to deter- 
mine the discriminative abilities of the three kinds of functional 
connectivity. Finally, we analyzed the abnormalities of non-linear 
and linear functional connections in schizophrenia and deter- 
mined the spatial distribution of the brain regions related to 
symptoms of schizophrenia. 

CLASSIFICATION 

We examined the whole brain resting-sate functional connectivity 
of the schizophrenic patients and the healthy controls using PCC, 
MIC, and eMIC, respectively. The 6670 mean functional connec- 
tions of the controls (Figure 1, the first row) and the patients 
(Figure 1, the mid row) according to the PCC, MIC, and eMIC. 
In the third row of Figure 1, the first two panels were presented 
the relationship between the MIC and PCC, and between the 
eMIC and PCC, where each red star represented one of the 6670 
mean functional connections over the patients, each blue cross 
represented one of the 6670 mean functional connections over 
the controls. The last panel in the third row of Figure 1 showed 
the mean and deviation of the 6670 mean function connections 
of the patients (red bar) and the controls (blue bar) which were 
depicted in the first two panels of the third row. Linear SVM was 
employed to differentiate the patients with schizophrenia from 
healthy controls using the whole brain functional connectivity 
extracted by the PCC, MIC, and eMIC. 

To evaluate the discriminative power of the PCC, MIC, and 
eMIC, we conducted multivariate classification. When the num- 
bers of connections used for classification changed, the classifica- 
tion accuracies changed accordingly. The subjects were classified 
by using the first connections, and 6670 accuracy rates were 
obtained. With this approach, the eMIC exhibited excellent per- 
formance and the highest classification accuracy rate among the 
three methods. When the number of features for classification 
is fixed, the obtained accuracy cannot fully reflect the discrim- 
inative power of the whole feature space. Here, we present the 
classification results using a wide range (50, 100, 150, . . ., 1000) 
of the first selected connections (Figure 2). Clearly, the classifica- 
tion accuracies (GR) of the eMIC were consistently higher than 
the accuracies of the other two methods across the full feature 
space. 

CHANGES OF FUNCTIONAL CONNECTIVITY 

From Figure 2, we can see that when no more than 150 
features were involved in the classification, the PCC, MIC, 
and eMIC gave the best performance. Then, if more fea- 
tures were involved, the classification accuracy of the three 
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FIGURE 1 | Comparison of mean functional connectivity of the PCC, 
MIC, and eMIC. The top row showed the mean functional connectivity 
maps of the healthy controls, the mid row showed the mean 
functional connectivity maps of the schizophrenic patients, the bottom 
row (the first two) showed the 6670 mean functional connections of 



the patients (red star) and that of the healthy controls (blue cross) 
according to MIC, and eMIC against that of the PCC, and (the last 
one) the comparison of the mean and the deviation of all the 6670 
connections according to PCC, MIC, and eMIC between the patients 
(red bar) and the controls (blue bar). 




number of Features 



FIGURE 2 | Classification accuracy rates relative to the number of selected classification accuracy. GR (the rate of all the subjects correctly classified), SS 
connections extracted by the PCC, MIC, and eMIC. The x-axis represents the (the rate of the patients correctly classified), and SC (the rate of the controls 
number of connections involved in classification; the y-axis represents the correctly classified) are all plotted. Colors represent the types of connection. 



methods decreased and the classification accuracy tend to be 
stable when the number over 200. Therefore, consensus func- 
tional connections of the first 200 features involved in the 
classification corresponding to the PCC, MIC, and eMIC were 
evaluated. Additionally, the number of features was identified 
in accordance with previous study (Dosenbach et al, 2010). 
Then, 113 consensus functional connections were obtained 



from each iteration of the leave-one-out cross-validation for 
the MIC and PCC, and 122 consensus functional connec- 
tions were identified for the eMIC. The consensus functional 
connections from the PCC, MIC, and eMIC were projected 
to a surface rendering of a human brain that was visual- 
ized with BrainNetViewer (http://www.nitrc.org/projects/bnv/) 
(see Figure 3). Further, we added a probabilistic atlas of the 



Frontiers in Human Neuroscience 



www.frontiersin.org 



October 2013 | Volume 7 | Article 702 | 5 



Su et al. 



Nonlinear functional connectivity in schizophrenia 




FIGURE 3 | Region weights and strength of the connections constructed 
with PCC, MIC, and eMIC. The connections are displayed in a surface 
rendering of a human brain. The thicknesses of the consensus connections 
during leave-one-out cross-validation are scaled by their strength (normalized 
mean of the rank orders in all iterations during the leave-one-out 
cross-validation). Connections whose strength were increased in the patients 
relative to the controls are shown in orange, and connections whose strength 



were decreased in patients are shown in light green. The ROIs related to the 
selected consensus connections are also scaled by their weights (sum of the 
weights of all connections to and from that ROD and displayed. The ROIs are 
color-coded according to the functional networks (CON, white; DMN, green; 
cerebellum, red; visual network, yellow; sensorimotor network, cyan; 
frontal-parietal network, rose; and other, black). The number in this figure for 
ROIs is shown in Table Al. 



cerebellum to the ICBM152 cerebrum template that was released 
with the software. 

Primary consideration was given to changes in the strength of 
the connections. The consensus functional connections according 
to PCC and MIC were all decreased in patients with schizophre- 
nia compared to healthy controls. It is generally accepted that 
the strength of functional connectivity is decreased in patients 
with schizophrenia compared with healthy controls. The PCC- 
and MIC-based functional connections found in our study are 
consistent with this common view. In contrast, the strength of 
the eMIC connections were elevated in patients with schizophre- 
nia. In the Discussion section, we provide a detailed explanation. 
Additionally, in terms of anatomical location, both the con- 
sensus MIC connections and the consensus eMIC connections 
were largely consistent with the consensus PCC connections. The 
MIC and the PCC shared 56 common connections, the eMIC 
and the PCC shared 53 common connections, and the common 
connections make up approximately half of the total consensus 



connections. Significant difference of eMIC to PCC is related to 
the functional connections between the cerebellum and the tem- 
poral cortex located in the default model network (DMN) that 
were found by eMIC but were not detected by PCC. Although 
there were no PCC functional connections between the cerebel- 
lum and the temporal cortex, the cerebellum and the temporal 
cortex were both connected with the parietal cortex (Figure 4). 

DISTRIBUTION OF FUNCTIONAL CONNECTIVITY 

In addition, we investigated the distribution of brain regions that 
are related to the consensus connections produced by the PCC, 
MIC, and eMIC. Overall, the distributions of ROIs identified by 
the PCC, MIC, and eMIC were similar, although the consen- 
sus connections were not completely identical across the three 
methods (Figures 3, 4). 

To facilitate analysis, these ROIs can be categorized as fol- 
lows: (i) the cingulo-opercular network (CON), including several 
regions in the anterior prefrontal cortex, the inferior parietal 
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cortex, the basal ganglion, the dorsal anterior cingulate cortex 
(dACC), the insular, the thalamus, and the cerebellum (which 
will be discussed later separately due to its importance), which 
is a newly defined cognitive network with great significance for 
schizophrenia (Dosenbach et al, 2007; Tu et al., 2012); (ii) the 
DMN, including structures of the hippocampus, the posterior 
cingulate cortex, the medial prefrontal cortex (mPFC), and the 
bilateral inferior parietal cortex, which is believed to play an 
important role in the pathogenesis of schizophrenia (Raichle 
et al, 2001; Fransson, 2005; Whitfield-Gabrieli et al, 2009); 
(iii) the cerebellum network, which can be seen as part of the 
CON; (iv) the visual network comprising the primary visual cor- 
tex, the extra-striate visual areas and the lingual gyrus, fusiform 
gyrus, and calcarine gyrus, which is involved in visual processing 
(Beckmann and Smith, 2005; Damoiseaux et al, 2006; van den 
Heuvel et al, 2008; van den Heuvel and Hulshoff Pol, 2010a,b); 
(v) the sensorimotor network, including the primary sensory cor- 
tex, the primary motor cortex the supplementary motor cortex 
(Biswal et al, 1995; van den Heuvel and Hulshoff Pol, 2010a,b); 
and (vi) the frontal-parietal network, consisting of the superior 
parietal and the superior frontal cortex, which is involved in atten- 
tion processing (Dosenbach et al., 2010; Beckmann et al., 2005; 
De Luca et al, 2006) (Figure 4). 

For the three kinds of functional connectivity, the ROIs with 
the heaviest weight were distributed primarily in the DMN, cere- 
bellum, visual network, and frontal-parietal network. Specifically, 
the cerebellum was the most important network according to 



eMIC, while the visual cortex was the most weighted network 
according to PCC and MIC. 

DISCUSSION 

In this study, we introduced a novel measure called eMIC for 
estimating the non-linear functional connectivity underlying 
schizophrenia and applied this estimation of functional connec- 
tivity to distinguish schizophrenic patients from healthy controls. 
Then, we found that strength of the non-linear functional con- 
nectivity increased in patients with schizophrenia which was 
opposed to that observed for the traditional method, which can 
be attributed to the compensatory mechanisms in the human 
brain. Furthermore, the non-linear and linear functional con- 
nectivity presented similar but not completely the same spatial 
distribution. 

ANALYSIS OF CLASSIFICATION 

The results of the classifications produced by the PCC, MIC, and 
eMIC were obtained using the same procedure and the same 
classifier parameters. The only factor affecting the classification 
accuracy is the measure of the functional connectivity. Using sup- 
port vector classification, we compared the discriminative powers 
of the three kinds of functional connectivity produced by PCC, 
MIC, and eMIC. The eMIC produced consistently higher clas- 
sification accuracies than the other two methods across the full 
connection space (Figure 2). Furthermore, when all 6670 features 
were used for classification, we obtained the classification results 
displayed in Table 2. The classification accuracy of the eMIC 
remained higher than those of the other methods. In short, the 
functional connectivity produced by the eMIC had the greatest 
discriminative ability among all three methods. 

Now, we pay more attention to the reason for the better classi- 
fication accuracy of eMIC. The MIC maximizes the association 
between two time series, whereas the PCC captures only lin- 
ear function. Then, the eMIC may capture the subtle non-linear 
neuronal synchronization in human brain, which will improve 
the performance of eMIC. From another point of view, the 
eMIC combines the discriminative information of both the PCC 
and MIC, which reflects a feature-level information fusion that 
increased classification accuracy. An Mi-based study reported 
that the linear functions accounted for most of the associa- 
tions between the fMRI time series and that non-linear functions 
only mined 5% more information in fMRI time series (Hlinka 
et al, 2011); this could be explain why the eMIC gave limited 
improvement to PCC (Figure 2). Conservatively, this result con- 
firmed that the non-linear associations have discriminative power 
which should not be overlooked. We believe that these non-linear 
connections add new discriminative information to the linear 
connections which may increase the classification accuracy. 



Table 2 | Classification results when all 6670 features were involved. 

Method (%) SC (%) SS (%) GR (%) 

MIC 71.9 81.2 76.6 

PCC 81.2 81.2 81.2 

eMIC 81.2 84.4 82.8 




FIGURE 4 | Distribution of selected consensus connections 
constructed by the PCC, MIC, and eMIC and region weights of related 
ROIs demonstrated in a circle graph. The names of the ROIs are 
color-coded as shown in Figure 3. The green lines represent connections 
constructed by the PCC, blue lines represent connections constructed by 
the MIC, and orange lines represent connections constructed by the eMIC. 
This figure was plotted using the MATLAB toolbox called PlotPie which was 
developed by our study group and will be released in the near future. 
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CHANGES OF FUNCTIONAL CONNECTIVITY IN SCHIZOPHRENIA 

Here, we give a detailed analysis of the increased strength of the 
selected eMIC functional connections. The decreased strength 
of the PCC functional connections in schizophrenia have been 
noted in an overwhelming majority of studies (Camchong et al., 
2011; Pettersson-Yeo et al, 2011; Repovs et al, 2011; Fornito 
et al., 2012), which was also confirmed by our previous study 
(Shen et al., 2010). However, the eMIC-based functional connec- 
tions demonstrated increased strength in schizophrenic patients 
compared to healthy controls. As MIC captures various associ- 
ations but the majority of associations were linear which can 
be estimated by PCC (Hlinka et al., 2011), decrease of the PCC 
strength leads to decrease of the MIC strength. The human 
brain is an organ with great plasticity and adaptability. Thus, the 
compensatory mechanism of the human brain, that's non-linear 
correlation will strengthened to reconcile the influence of the cor- 
rupt of linear correlation, can result in increases of non-linear 
functional connectivity. Although the connection selecting pro- 
cedures of the three methods were implemented separately, the 
consensus functional connections by the eMIC and the PCC still 
had a relatively high identity (approximately 50%), which sup- 
port the notion that the compensatory mechanism of the human 
brain results in the increase strength of the eMIC in patients with 
schizophrenia. 

If the strength of MIC did not change between the patients and 
the controls, while the PCC decreased, the increase of the strength 
of eMIC will be purely caused by the changes of PCC. In fact, the 
connection strength of MIC decreased in the patients (Figure 1) 
and gave classification accuracy over 80% (Figure 2). Thus, the 
increase of eMIC in schizophrenic patients was not simply caused 
by the decrease of PCC but by the decrease of both PCC and 
MIC. We believe that the eMIC captures subtle changes of func- 
tional connectivity and adds new discriminative information to 
the classification. 

NETWORK ANALYSIS OF NON-LINEAR FUNCTIONAL CONNECTIVITY 

As the eMIC provided useful discriminative information to the 
linear functional connectivity, aberrant non-linear and linear 
functional connections were both categorized into six networks 
for explaining the symptom of schizophrenia. Therein, connec- 
tions between the DMN and the cerebellum, the DMN and the 
visual network, the DMN and the frontal-parietal network, and 
the cerebellum and the frontal-parietal network demonstrated the 
greatest discriminative power for both the linear and non-linear 
functional connectivity. 

The DMN, the CON (cerebellum), visual network, and 
frontal-parietal network which are related to specific brain func- 
tions were all weighted by the linear and non-linear measures. The 
DMN was frequently reported in previous studies and weighted in 
our study. The three methods all identified the important nodes 
in the DMN such as the mPFC, inferior parietal cortex, Para 
hippocampal gyrus, and middle temporal cortex. This network 
is generally accepted as an important network associated with 
schizophrenia (Bluhm et al., 2007; Garrity et al., 2007; Whitfield- 
Gabrieli et al., 2009). The CON is believed to support the "task 
model" in opposition to the "default model" (Dosenbach et al., 
2006). In our study, the cerebellum_6_R which can be seen as part 



of the CON is important for all the three functional connectiv- 
ity, particularly for the new non-linear measure. Additionally, the 
PCC, MIC, and eMIC all gave the visual network heavy weight. 
The frontal-parietal network was also identified by the all the 
three measures, which was reported as important brain regions 
in previous studies (Honey et al, 2005; Lynall et al., 2010). 

The DMN, the CON (including cerebellum), visual network, 
and frontal-parietal network obtained slightly different region 
weight according to the PCC, MIC, and eMIC, respectively. 
Functional connections between the DMN and the cerebellum 
were identified by the eMIC but not by the PCC. If the cere- 
bellum is viewed as part of the CON, interactions between 
the DMN and the CON can be established by the eMIC. The 
wide spread changes in the CON and DMN may interact with 
each other, and contribute to the functional basis of schizophre- 
nia. In the connections produced by the PCC and MIC, the 
right inferior occipital cortex in the visual network exhibited 
the greatest region weights. However, the cerebellum_6_R exhib- 
ited the greatest region weight for the eMIC. This result implies 
that the cerebellum may play a role in non-linear interac- 
tion between different brain regions. Connections between the 
visual network and regions of the frontal cortex are believed 
to be involved in visual perception which may contribute to 
the aberrance of visual perception in schizophrenia (Harvey 
et al, 2011; Calderone et al., 2013). The cerebellum and the 
DMN are linked to the frontal-parietal network by eMIC and 
MIC but not PCC. Combined with the DMN, this network is 
thought to be closely related to attention tasks (Bush and Shin, 
2006; Gao and Lin, 2012), which may explain the occurrence 
of attention impairment and its relationship to the DMN in 
schizophrenia. 

In conclusion, from the functional network perspective, the 
distribution of ROIs with the greatest weight according to the lin- 
ear and non-linear connections was similar but not completely 
the same, and the non-linear connections shed new light on 
interpretation of the schizophrenic symptoms. 

LIMITATIONS 

Although the new functional connectivity constructed by using 
the eMIC method exhibits better performance in the classification 
and has explored new information about schizophrenia, there are 
several limitation in this study. First, our sample size was small. 
Generalizations of the proposed methods need to be tested with 
large data sets. Second, the number of slices in the fMRI image 
used in this study was 20, which was relatively fewer and may 
not sufficiently to capture the details of the abnormalities in the 
patients. Third, network analysis of the human brain is a trend 
in the literature which include ROI localization and connectivity 
estimation. Our study is confined to connectivity analyses that do 
not include definition of ROIs. 

CONCLUSIONS 

In this study, we introduced a novel non-linear functional con- 
nectivity for schizophrenia study. The classification results show 
that the non-linear functional connectivity has an equal if not 
better discriminative ability than existing linear measures in 
schizophrenia identification. This result suggests that non-linear 
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functional connectivity should be taken into account in 
research on schizophrenia and other psychiatric disorders. 
Furthermore, we found that the non-linear functional con- 
nectivity which was strengthened in the patients has a 
similar distribution with its linear counterpart. This new 
finding indicates the presence of compensatory mechanism 
between linear and non-linear associations and non-linear 
functional network abnormalities underlying the symptom of 
schizophrenia. 
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APPENDIX 

Table A1 | Names of the ROIs used in this study. 

1. Amygdala_L 

2. Amygdala_R 

3. Angular_L 

4. Angular_R 

5. Calcarine_L 

6. Calcarine_R 

7. Caudate_L 

8. Caudate_R 

9. Cerebelum_3_L 

10. Cerebelum_3_R 

11. Cerebelum_4_5_L 

12. Cerebelum_4_5_R 

13. Cerebelum_6_L 

14. Cerebelum_6_R 

15. Cerebelum_7b_L 

16. Cerebelum_7b_R 

17. Cerebelum_8_L 

18. Cerebelum_8_R 

19. Cerebelum_9_L 

20. Cerebelum_9_R 

21. Cerebelum_10_L 

22. Cerebelum_10_R 

23. Cerebelum_Crus1_L 

24. Cerebelum_Crus1_R 

25. Cerebelum_Crus2_L 

26. Cerebelum_Crus2_R 

27. Cingulum_Ant_L 

28. Cingulum_Ant_R 

29. Cingulum_Mid_L 

30. Cingulum_Mid_R 

31. Cingulum_Post_L 

32. Cingulum_Post_R 

33. Cuneus_L 

34. Cuneus_R 

35. Frontal_lnf_Oper_L 

36. Frontal_lnf_Oper_R 

37. Frontal_lnf_Orb_L 

38. Frontal_lnf_Orb_R 

39. Frontal_lnf_Tri_L 

40. Frontal_lnf_Tri_R 

41. Frontal_Med_Orb_L 

42. Frontal_Med_Orb_R 

43. Frontal_Mid_L 

44. Frontal_Mid_Orb_L 

45. Frontal_Mid_Orb_R 

46. Frontal_Mid_R 

47. Frontal_Sup_L 

48. Frontal_Sup_Media_L 

49. Frontal_Sup_Media_R 

50. Frontal_Sup_Orb_L 

51. Frontal_Sup_Orb_R 

52. Frontal_Sup_R 

53. Fusiform_L 

54. Fusiform_R 

55. HeschLL 

56. HeschLR 

57. Hippocampus_L 

(Continued) 
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58. Hippocampus_R 

59. lnsula_L 

60. lnsula_R 

61. Lingual_L 

62. LinguaLR 

63. Occipital_lnf_L 

64. Occipital_lnf_R 

65. Occipital_Mid_L 

66. Occipital_Mid_R 

67. Occipital_Sup_L 

68. Occipital_Sup_R 

69. Olfactory_L 

70. Olfactory_R 

71. Pallidum_L 

72. Pallidum_R 

73. Paracentral_Lobule_L 

74. Paracentral_Lobule_R 

75. ParaHippocampal_L 

76. ParaHippocampal_R 

77. Parietal_lnf_L 

78. Parietal_lnf_R 

79. Parietal_Sup_L 

80. Parietal_Sup_R 

81. Postcentral_L 

82. PostcentraLR 

83. PrecentraLL 

84. PrecentraLR 

85. Precuneus_L 

86. Precuneus_R 

87. Putamen_L 

88. Putamen_R 

89. Rectus_L 

90. Rectus_R 

91. Rolandic_Oper_L 

92. Rolandic_Oper_R 

93. Supp_Motor_Area_L 

94. Supp_Motor_Area_R 

95. SupraMarginaLL 

96. SupraMarginaLR 

97. Temporal_lnf_L 

98. Temporal_lnf_R 

99. Temporal_Mid_L 

100. Temporal_Mid_R 

101. Temporal_Pole_Mid_L 

102. Temporal_Pole_Mid_R 

103. Temporal_Pole_Sup_L 

104. Temporal_Pole_Sup_R 

105. Temporal_Sup_L 

106. Temporal_Sup_R 

107. Thalamus_L 

108. Thalamus_R 

109. Vermis_1_2 

110. Vermis_3 

111. Vermis_4_5 

112. Vermis_6 

113. Vermis_7 

114. Vermis_8 

115. Vermis_9 

116. VermislO 
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